Probability Distribution Function of the Order Parameter: Mixing Fields and 

Universality 
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We briefly review the use of the order parameter probability distribution function as a useful 
tool to obtain the critical properties of statistical mechanical models using computer Monte Carlo 
simulations. Some simple discrete spin magnetic systems on a lattice, such as Ising, general spin-S 1 
Blume-Capel and Baxter- Wu, Q-state Potts, among other models, will be considered as examples. 
The importance and the necessity of the role of mixing fields in asymmetric magnetic models will 
be discussed in more detail, as well as the corresponding distributions of the extensive conjugate 
variables. 

PACS numbers: 75.10Hk,64.60.Ht,75.40.Mg 



I. INTRODUCTION 

According to the words of Sir Humphry Davy (1778- 
1829) [l|, "Nothing tends to the advancement of knowl- 
edge as the application of a new instrument. " As a mat- 
ter of fact, several examples can be traced back to 
levers, thermometers, telescopes and, perhaps nowadays, 
to computers. In particular, the use of the thermometer, 
developed during the 16th and 17th centuries, amongst 
other fantastic applications, also marked the dawn of the 
thermal physics and phase transitions era (not count- 
ing, of course, the discovery of fire itself). Apart from 
its rapid technical application in medicine, it is believed 
that the thermometer in thermal physics can be rightly 
compared to the use of the telescope in astronomy. It 
also made possible the discovery of the specific heat and 
latent heat by Joseph Black (1728-1799) in Scotland (see, 
for instance, references and @|). On April 1762, Black 
coined the term latent heat to refer to the amount of en- 
ergy required to change the phases solid <-> liquid or liq- 
uid <H> vapor without changing their temperature. This is 
nowadays known as a first-order, or discontinuous, phase 
transition. 

In 1822, Charles Cagniard de la Tour (1777-1859) 
in France discovered critical phenomena 0- Although 
Cagniard actually saw the liquid and vapor fluid phases 
becoming equal, the term critical was only coined some 
years later, in 1859, by Thomas Andrews (1813-1885) 
through the observation of what he called critical opales- 
cence, first observed in carbon dioxide (carbon dioxide 
was discovered by the same Joseph Black almost one 
century before). At this critical point there is neither 
distinction between the liquid and vapor phases nor la- 



tent heat, and such a transition is nowadays known as 
a second-order, or continuous, phase transition (for a re- 
cent review on phase transitions see Kadanoff [J], and 
also Domb @). 

Together with the advent of these completely new phe- 
nomena, a question was raised: the determination of the 
critical point for different substances. According to the 
Gibbs phase rule, a single component system can present 
only isolated critical points, as well as only single triple 
points [6J . Although the critical temperature T c and the 
critical pressure p c are rather easily measured, the mea- 
surement of the critical density p c is not so straightfor- 
ward and turns out to be a much more difficult task. 
This can be seen in Figure [1] from the experimental re- 
sults of three different compounds. As an example, in this 
Figure one has the original data reported by Louis Paul 
Cailletet (1832-1913) and Emilc Ovide Joseph Mathias 
(1861-1942), where the lack of data to measure p c is 
apparent Q- It was actually Cailletet and Mathias who 
proposed what was later (and still is) known as the law 
of the rectilinear diameter 



Pd = t;(pl + Pv) = Po + AT, 



(1) 
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where pd is the diameter density, p^ and pv are the liq- 
uid and vapor densities, respectively, and p and A are 
constants which depend on the substance. The critical 
density is then obtained from the measurement of T c as 
Pc = Po + AT C . The reader interested in the further 
developments and details of the empirical law of the rec- 
tilinear diameter is directed to the recent review by Reif- 
Acherman [8|. 

Besides the study of the phase diagram of several sub- 
stances, in which triple and critical points have been de- 
termined, the phase transitions of some magnetic mate- 
rials were also well established by the end of the 19th 
century. In fact, in 1895, Pierre Curie (1859-1906) put 
the correspondence between magnets and fluids on a solid 
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FIG. 1: Original data of the proposed law of rectilinear di- 
ameter by Cailletet and Mathias for nitrous oxide, ethylene, 
and carbon dioxide (taken from reference [?()• In the abscissa 
one has the temperature in the Celsius scale and in the or- 
dinate the density in 10 -3 kg/m 3 . For each compound, the 
upper curve is the density of the liquid, and the lower curve 
the density of the vapor. 




FIG. 2: Sketch of the magnetization M as a function of tem- 
perature T, at zero external field, for a simple magnet. M is 
the saturation magnetization and T c is the critical tempera- 
ture. 



basis by taking the pressure p as the analogue of the ex- 
ternal magnetic field H , and the density p as the analogue 
of the magnetization M [5j. As a simple example, the 
sketch in Figure [5] for the magnetization as a function 
of temperature for a simple symmetric magnet clearly 
shows the similarity to the density versus temperature 
for a fluid in Figure [TJ By symmetric here one means 
that, at zero external field, the free energy is the same 
whether the magnetization points in one direction or the 
other. This symmetry for the magnet, however, implies 



that the equivalent of the diameter for the magnetiza- 
tion is temperature independent, and any theory applied 
close to the critical point should be able to describe such 
diverse behavior. 

An astonishing discovery came quite later regarding 
the universal behavior of different substances close to a 
continuous phase transition. Measurable thermodynamic 
quantities, such as the order parameter m (which can be 
taken as the magnetization M for the magnet and the 
difference pl — Pv for the fluid), specific heat at constant 
volume cy, magnetic susceptibility \ ( or the compress- 
ibility in the fluid case), present a power law behavior of 
the form 
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where means asymptotically equal to. In the above 
equations, m , C a , Xo, and mi are constants, /3, a, 7, 
and 5 are the corresponding thermal critical exponents, 
and t = T — T c (for questions of clarity, and without 
loss of generality, in what follows we will use this defini- 
tion of t instead of the more usual one t = (T — T c )/T c ). 
More (spatial) exponents can be defined for the correla- 
tion function (exponent 77), correlation length (exponent 
v), and also for the dynamical critical slowing down (dy- 
namical exponent z). However, in the present work, for 
simplicity, we are going to treat only the thermal expo- 
nents in Eqs. @. What happens is that, regardless of 
the microscopic structure of the interactions among the 
constituents of the compounds, the critical exponents are 
the same for systems belonging to the same universality 
class. A universality class, according to Kadanoff, can be 
defined primarily by the spatial dimension of the lattice, 
the symmetry of the order parameter and the range of in- 
teractions 0, • One has also to consider the number of 
phases becoming equal at the continuous transition. In 
addition, the scaling laws also provide relations among 
critical exponents in such a way that knowing just two 
of them allows one to determine all the remainders. For 
example, just to cite a couple of them due to the works 
of Rushbrooke and Widom, we have (ill 



a + 2/3 + 7 = 2; 
136 = + 7. 



(Rushbrooke) 
(Widom) 



(3) 



A theoretical approach which could give a satisfactory 
picture of the scaling laws was proposed by Widom and 
some other workers (see, for instance, reference Q and 
references therein). The Widom scaling hypothesis just 
states that the singular part of the free energy is a gener- 
alized homogeneous function of its variables. From this 
hypothesis all the scaling laws follow, but the general- 
ized homogeneous function does not provide the values 
of the critical exponents. The renormalization group 
(RG) approach of Wilson [l2|, could give a complete 
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microscopic account for the Widom hypothesis and the 
Kadanoff universality ideas, as well as turned out to be 
a quite powerful way of getting numerical values for the 
critical exponents. The RG can indeed be considered a 
cornerstone in the study of critical phenomena. 

It should also be stressed that not only are the critical 
exponents universal but so are the scaling functions of the 
corresponding thermodynamic variables. For example, 
the universal character of the order-parameter probabil- 
ity distribution function (PDF) has been independently 
introduced by Bruce [l4| and Binder in 1981 [Hj. Its 
use in describing critical properties of models in statisti- 
cal mechanics has been shown to be quite valuable and 
has also been extended in the study of pure and disor- 
dered magnetic systems, e.g. Lennard- Jones fluids, crit- 
ical point in the unified theory of weak and electromag- 
netic interactions and quantum chromodynamics, non- 
equilibrium phase transitions, among others. Due to the 
universal behavior of the order-parameter PDF, it be- 
comes a fingerprint of the criticality of each class of sys- 
tems. If a match with this distribution is achieved, not 
only are the critical exponents in the same universality 
class, but also non-universal behavior (such as the deter- 
mination of the corresponding critical temperature) can 
be readily extracted. 

Without a doubt, one of the best ways to compute 
the order-parameter PDF presented above is the use of 
Monte Carlo (MC) computer methods[16j. As was said 
previously in this Introduction, the application of the 
new instrument called computer has in some sense rev- 
olutionized not only Physics but the sciences in general. 
For more than half a century, computer simulations have 
been shown to be useful, and even crucial, to the under- 
standing of most features and important subjects in mod- 
ern science. Indeed, in what concerns the present subject, 
MC simulations have been proven to be extremely valu- 
able for obtaining information about phase transitions, 
regarding their universal and non-universal aspects, such 
as critical exponents, critical amplitudes, critical tem- 
perature, and even correction-to-scaling properties. The 
essential idea behind an MC procedure is to simulate a fi- 
nite real system by generating a sample of states in which 
it can be found. However, for a very large system, only 
a small fraction of the total number of possible states 
is sampled. This leads, of course, to the inconvenience 
of having statistical errors associated with the physical 
quantities that are evaluated during the simulation. Nev- 
ertheless, with the fast increase of computer power in 
recent decades and the development of more efficient al- 
gorithms, high-precision estimates can be achieved with 
sufficiently large runs. Furthermore, as computer sim- 
ulations consider only finite systems, one must have a 
mechanism to extract the critical behavior (in the ther- 
modynamic limit) from the size dependence of the mea- 
sured physical quantities. In other words, it should be 
possible to get the non-analytical behavior of the infinite 
system at criticality by just studying the corresponding 
analytical behavior of finite systems. It turns out that 



such a procedure, in some cases, works quite well even for 
systems not too large. According to the finite-size scal- 
ing theory, proposed by Fisher in 1971, since near the 
transition the correlation length is limited by the system 
size L, one can get a set of relations similar to Eqs. ([2]) 
which includes the dependence on L as well [13, [HI]. For 
instance, at criticality (T = T C ) the extension of Eqs. © 
for the finite-size effects of the systems can be written as 

tow m'eL-P/", 

X« X' ^'\ (4) 

where m' , C' a , and x'o are constants. For any MC simu- 
lations, the finite-size scaling must be taken into account 
to avoid misleading results. Furthermore, one has also 
to spend some time in evaluating the errors coming from 
the simulations and the time correlations between config- 
urations due to the pseudo-random number generators. 

We will herein discuss the usefulness of the probabil- 
ity distribution function of the order parameter, or any 
other conjugate extensive variable, in obtaining critical, 
and mainly multicritical, behavior of magnetic physical 
systems. Emphasis will be given on Monte Carlo simula- 
tions to compute the order-parameter distribution func- 
tion, or any other needed distribution. We will discuss 
symmetric and (more ubiquitous) asymmetric cases. In 
order to get insight of the role of asymmetry in some 
models we will start with a brief revision of the mixing 
fields for fluid systems in the next section. We will closely 
follow the revised scaling by Rehr and Mermim [l!| and 
the corresponding field mixing procedure introduced by 
Bruce and Wilding [2(1 Hl[ , since for the present systems 
it is not necessary to consider the general approach of 
the complete scaling proposed by Kim and Fisher in the 
study of highly asymmetric fluids [22T - [24| . The order pa- 
rameter PDF in the symmetric case with its finite-size 
scaling relation will be derived in section III. In section 
IV we will apply it to the simplest Ising and Baxter- Wu 
models. A generalization of the finite-size scaling relation 
taking into account the presence of asymmetric fields will 
be discussed in section V. Examples of discrete spin sys- 
tems as the Blume-Capel and spin-1 Baxter- Wu models, 
together with the Lennard- Jones fluid, will be discussed 
in section VI. In the last section, the use of the proba- 
bility distribution to some other related models will be 
briefly presented together with some final remarks. 

II. FREE ENERGY SCALING FUNCTION 

For any symmetric system it is well known that the 
free energy g(T, A), in terms of the temperature T and 
the field A, can be written as 

g(T, A) = 9o (T, A) + \tf- a f ± (A/\tf s ), (5) 

where g (T, A) is an analytic function, t = (T — T c ) and 
the scaling function / + holds for T > T c and /_ for 
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T < T c , with the exponents defined in Eqs. ([2]). The 
symmetry here means that the critical point is located 
at t = Q and A = 0. This is the case for the Ising model 
where we have the correspondence A — H , H being the 
external magnetic field, as depicted in Figure [3] 







FIG. 3: Sketch of the symmetric phase diagram in the tem- 
perature field plane. The full circle is the critical point. First- 
order phase transitions take place for T < T c and A = 0. This 
is the equivalent field phase diagram of Figure [2] 

However, it is more common to have a first-order tran- 
sition line where the field A is a function of the tempera- 
ture T, namely A = A(T). This is the case, for instance, 
in the Blume-Capel model where A plays the role of the 
crystal field (or single ion anisotropy), and in fluid sys- 
tems where A can be the pressure or chemical potential. 
A typical phase diagram is sketched in Figure 0] and a 
possible generalization of the scaling function ([5]) can be 
written as 

g(T, A) = 9o (T, A) + \t\ 2 ~ a f ± [(A - A(t))/\tf s ], (6) 



A(T) 



A-A r 




T-Xc 



FIG. 4: Sketch of the first-order transition line A(T) in the 
field plane. (T c , A c ) is a critical or, in general, a multicritical 
point when one has an additional second-order transition line. 
t and rj are the field mixing directions. All points (T, A) along 
the r direction (such as point A) have rj — while points like 
B have r = 0. 



([5]) should be compatible with the power law behavior 
of the order parameter and compressibility given in Eqs. 
([2]) as well as with the rectilinear diameter law in (Q~]). 

Let us for the moment consider g = p and A = p in 
order to treat in more detail the fluid system according 
to ([B]). As discussed in reference [llj], for T < T c the 
function f—(z) is an even function of the variable z with 
the liquid phase holding for z — >• + and the vapor phase 
for z — > 0~ . We then have 



(||) =/<-+ 'I ''/-('-)• 



(7) 



The above scaling function, as well as its extension 
by including mixing fields, has been extensively stud- 
ied by Rehr and Mermin [19fl in the context of the fluid 
critical point. In reference [ijj], instead of the free en- 
ergy, the pressure, as a function of temperature and 
chemical potential, has been considered. As the Gibbs 
free energy per particle of a fluid is the chemical po- 
tential (i written in terms of the intensive temperature 
and pressure variables, i.e. p, = p(T,p), this equation 
can be inverted in order to get p = p(T,p). In doing 

this one gets equations of state of the form ^ = 

and p = i = ^f^^ , where the entropy per particle 

is s = S/N, v = V/N is the volume per particle, and 
p is the (number) density. In addition, the isothermal 

compressibility kt — — y I €~ ) can also be given by 




This implies that the scaling function 



^=(10 =1*™*), (8) 

where we have used Eqs. ([3]) in order to get 2 — a — (35 = (3 
and 2 — a — 2(35 — —7. It is clear from Eq. © that the 
compressibility indeed satisfies the corresponding power 
law behavior expressed in Eqs. (f2J). For the density Eq. 
© gives 

PL=Pc + \tffL(0 + ), Pv = p c +\tff_(o-), (9) 

and as f'_(0+) = -/i(0") one has 

p L -p v = 2f'_(0+)\tf, (10) 

meaning that the fluid order parameter pl~ Pv also pos- 
sesses its expected scaling behavior. However, from Eqs. 
© one cannot obtain the law of rectilinear diameter in 
(fTJ) unless A = 0. This is alright for the symmetric case 
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depicted in Figure [2] where the corresponding diameter 
of the magnetization is temperature independent. Al- 
though corrections to scaling in Eq. (j6|) can provide the 
expected behavior for the diameter in asymmetric sys- 
tems, Rehr and Mermin [l9j have proposed a field mix- 
ing procedure which is not so drastic as the corrections to 
scaling. Instead of working with the natural fields T — T c 
and A — A c , as illustrated by the straight dashed lines 
in Figure [4j one can define more appropriate fields in the 
following way. 

Consider the first-order transition line in the T x A 
plane in such a way that it is given by a function A(T). 
Close to the multicritical point one has 

A(T) = A c - r(T - T c ) + (11) 

where r = tan^f (see Figure SJ. If we neglect the higher- 
order terms in Eq. fTTj) we have 

(A - Ac) + r(T - T c ) = 0. (12) 

The equation above gives a straight line tangent to the 
A(T) curve at (T c , A c ), so we can define a field 

77 = (A-A c ) + r(T-T c ), (13) 

which may have any direction not parallel to the tangent 
line. Similarly, we have 

t = (T - T c ) - s(A - Ac), (14) 

where s = tan#. 7/ = defines the r axis while r = 
defines the 77 axis, as shown in Figure |4j When r = s = 
one gets the symmetric choice (r = T — T c , 77 = A — A c ). 
Thus, in these new fields Eq. ((HJ) takes the form 

^,A)=. go (r,77) + |r| 2 -"/ ± (77/|r|^). (15) 

It should be noted that in the language of a fluid 77 — > + 
means the liquid phase while 77 — > 0~ means the vapor 
phase. From Eq. (TT51) one gets 

PL-pv=2p-f'-(0 + )\tf, (16) 



p d = p c + (2-a)^/_(0+)|r| 1 -«. (17) 

This means that the scaling relation (fT"5"|) indeed yields 
the singularities observed in the order parameter and in 
the diameter. Regarding the diameter itself one sees that 
the temperature exponent is not 1. As a is small the dif- 
ference is not so large. Nevertheless, some experiments 
do show deviations from Eq. (fTJ) and quite a few dif- 
ferent empirical approaches have been proposed to deal 
with them (the interested reader should consult refer- 
ence [7|). Moreover, the singularities of the isothermal 
compressibility, specific heat, as well as the difference 

(ffr) L ~~ (tOv x which are omitted here for 



simplicity, are also reproduced by the scaling relation 
with the mixed fields in Eq. (|15l) . It should be said 
that the latter difference is not obtained from the sim- 
ple relation © either [l9j]. Additional contributions to 
Eq. (fT7|) are obtained when one considers the complete 
scaling procedure^. 

Thus, whenever one deals with models having asym- 
metries in the phase diagram, one has to work with the 
mixing fields given by Eqs. (|13[) and (fT4|) rather than 
(t, A). This procedure can also be understood as a way 
of minimizing the corrections to scaling, always present 
in experiments for not being close enough to the critical 
or multicritical point, as well as in simulations for not 
being able to simulate large enough lattices. 

In what follows we will discuss these ideas with the use 
of the order-parameter probability distribution function 
in getting not only the universality class of the system 
but also locate the multicritical point, which is really a 
difficult task in some cases. 



III. ORDER PARAMETER PROBABILITY 
DISTRIBUTION FUNCTION 

For an infinite system undergoing a continuous phase 
transition at a temperature T c one has the non-analytic 
behavior expressed by Eqs. @ for those thermodynamic 
quantities. On the other hand, for a finite system of 
linear dimension L, with a number of sites N = L d , 
where d is the dimension of the lattice, those singulari- 
ties are smeared out and an analytic behavior is observed 
instead. For rather small values of L exact results (parti- 
tion function) can be obtained. However, as L increases 
the number of states of the system gets drastically larger 
in such a way that exact results become impractical, in 
some cases, even for system sizes of around a dozen par- 
ticles. In this way, Monte Carlo (MC) simulations turn 
out to be one of the most efficient methods to measure 
the thermodynamical variables like the magnetization or 
order parameter m of a magnetic system, and also its 
probability distribution Pi(m,T). 

At the critical temperature and for very large system 
sizes (L >> 1) one expects the following scaling relation 

Mm 

P L (m) = bL^ v P*(bL^ v m), (18) 

where v is the correlation length critical exponent, b is a 
non-universal metric constant, and P* is a universal func- 
tion. P* is a fingerprint of the corresponding universality 
class. Note that in the above equation m is considered 
a continuous variable. The task of computing P* from 
Monte Carlo simulations and, in addition, its utility in 
obtaining either the criticality or the multicriticality of 
symmetric and asymmetric models, will be exemplified 
in the next sections. 
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IV. PROBABILITY DISTRIBUTION IN THE 
SYMMETRIC CASE: ISING AND BAXTER- WU 
MODELS 

Let us start with the symmetric Ising model which can 
be defined by the following Hamiltonian 



n 



H 



(19) 



1=1 



where the first sum is over nearest-neighbor pairs of spins 
and the last one is over all N = L d sites of a hyper 
cubic d-dimensional lattice with linear size L. J is the 
exchange interaction and H is the external magnetic field. 
The spins Si can take values —S, —S + 1, S — 1, S. The 
simplest case is the spin-1/2 model where S = 1/2. 

For a finite-size system of linear dimension L and a 
Monte Carlo run with Af Monte Carlo steps per spin (af- 
ter discarding some initial configurations for equilibrating 
the system) one has, for each generated configuration k, 
a measure of the energy at zero field P£ and the corre- 
sponding magnetization 



N 



(id) 



(20) 



i=i 



in such a way that the Ising Hamiltonian (|19[) can be 
written as Hl = —J El - HM^. One then has the fol- 
lowing mean values for El, Ml and Hl 

k—l k—1 

(H L ) = -J {E L ) - H (M L ) . (21) 



For this run one can also compute Pt(Mi), which is 
the probability of having the magnetization Mj in this 
particular simulation (in this case —SN < Mi < SN) 

„ , , number of occurrences of Mi 

Pl {Mi) = 1 ^ J - l - . (22) 

Note that the mean value of Ml is also given by 

(M L ) =J2MiP L {Mi), (23) 



as expected, and from Eg. (122)) one also has the normal- 
ization condition 



Y j Pl(M 1 ) = 1. 



(24) 



In the more general case one can obtain the joint prob- 
ability distribution Pl(P(, Mi). The distribution given 
by Eq.|J22J) is obtained from P L (M l ) = P L (E l ,M i ), 
where -S 2 Nc/2 < E t < S 2 Nc/2, with c the coordina- 
tion number of the lattice. 



The next step is now to compute Pl {m) from the mea- 
sured Pi(Mj). Eq. (l2~4"]) can be rewritten as 



E 



Pl(M) 
AM 



AM = 1, 



(25) 



where AM = 1 is the discrete magnetization step (if, as it 
is usual for the spin-1/2 Ising model, one defines Si = ±1 
one would have AM = 2). By defining rm = J^>, one 
has a discrete series of non-integer values — 1 < m j < 1 , 
so the above equation can be written as 

\- PlK) AM _ x - P L (mQ _ 

^ AM/SN SN ~ 2^ AM / SN ~ [ ° J 

rrii — — l mi — — l 

Now, for a large system the quantity drrii — AM is very 
small and the magnetization m, can be regarded as a 
continuous variable so the above sum can be replaced by 
the integral 



PL(m)dm, = 1, 



where 



P L (m) 



PUmi) 
AM/SN' 



(27) 



(28) 



which is the desired expression for the continuous prob- 
ability distribution. 

This probability distribution satisfies the scaling rela- 
tion HH) 



P L (m) = bL l5 / v P* L {bL p /"m) = bL^" P* L {m*), (29) 

where m* = bL^l u m. 

The normalization condition 

J P L (m)dm = 1 = J bL fi/y Pl{m*)dm = J P£(m*)rfm* 

(30) 

implies that P£(to*) is already normalized. 
On the other hand, the variance is given by 



m 2 PL(m)dm : 



(31) 



m Pl(m*)dm* 



b 2 L 



2B/u 



b 2 L 2/3/ V ■ 



m ^L\rn jam , 
(32) 

(33) 



Thus, also normalizing the variance a* = 1, one has 
bL^l v = 1/cr and m* = mja and one gets rid of the 
scaling exponent 



P£(m/(j) = aP L {m), 



(34) 



where m, Pl(to) and a are measured in an MC simu- 
lation. In this way we obtain the universal normalized 
function P* with unit variance. 
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A. One-dimensional Ising Model 

Exact results can be obtained for the one-dimensional 
Ising model [l4j]. On the other hand, analytic results can 
be obtained as well in higher dimensions by using some 
approximate approach [14] . In the present case, however, 
we will be mostly concerned with the order-parameter 
PDF computed from Monte Carlo simulations. Such a 
procedure has been shown to be very useful not only for 
the PDF itself but also in order to get the transition 
point. Details for the exact PDF for the one-dimensional 
model using the transfer matrix technique can be found 
in reference 1141 ■ 



B. Two-dimensional Ising Model 

The two-dimensional universality class of the PDF for 
the Ising model has been studied by Binder [l5} and Nico- 
laides and Bruce [13] ■ As a matter of illustration, Figure 
[5] shows the probability distribution of Pi(mj) according 
to Eq. (j2"2")h where the order parameter has been normal- 
ized by rrii = Mi/NS. In this simple example one has 
a rather small latice L — 32. The method used was the 
Metropolis algorithm for spin- 1/2, spin-1, and spin-3/2 
with just 6 x 10 6 MCS per spin. The inset shows the 
corresponding continuous limit of the probability distri- 
bution resulting from taking Eq. (|28[) . Different curves 
are obtained for different values of S. Finally, Figure 
IH1 depicts the two-dimensional Ising probability distribu- 
tion function of the order parameter given by Eq. Q34p. 
where one can clearly see the collapse of all the curves. 
For comparison, in Figure [5] one also has the data from 
the Wolff algorithm applied to the spin-1/2 Ising model 
for a lattice size L — 100 with 10 8 MCS. One can notice 
that, in this case, there is a very definite universal func- 
tion which is independent of the spin value, as expected. 
Larger lattices and longer MCS will provide quite similar 
results with smaller finite-size corrections and smoother 
curves. 

Figure [5] clearly shows the universal character of the 
probability distribution function P* , which could be com- 
puted by knowing, a priori, the transition temperature. 
In spite of that a different approach has been proposed to 
get not only the critical temperature but also the order- 
parameter PDF in cases where one does not know the 
exact location of the continuous transition 30] . One first 
generates the PDF PL (m,T ), where T a , although close 
to the actual transition temperature, can be either above 
or below T c . Simulations are then performed to obtain a 
matching of Pj,(m, T£) with PL o (m,T ) for L ^ L a at a 
temperature Tl ^ T Q . A usual finite-size scaling analy- 
sis is then used to get T c and the exponent v, which in 
turn allows one to get the desired universal function P*. 
This has been applied to the two-dimensional spin-1/2 
and spin-1 Ising models [30j . 

Knowledge of P* for this universality class can also 
be very useful in treating the diluted two-dimensional 
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m. 

i 

FIG. 5: Probability distribution function P(rrii), where 
rrii — Mi/NS, according to Eq. (|22[) for the two-dimensional 
Ising model with spin-1/2, 1, and 3/2. The data correspond 
to lattice size L = 32. The exact T c — 2.2692... has been 
taken for the spin-1/2, T c = 1.6935 for the spin-1 [2t| and 
T c = 3.2879 for the spin-3/2 The inset shows the corre- 
sponding distribution P(m) in the continuum limit from Eq. 
(|28[) for the same lattice and spin values. The error bars are 
smaller than the symbol sizes. 



Ising model. In this case, the Harris criterion is not ful- 
filed and there is not yet a consensus whether or not the 
critical exponents change with dilution. By employing 
MC simulations and using P* as a guide, it was possible 
to analise the corresponding order-parameter PDF [3lj . 
The results agree with the strong universality scenario, 
which predicts only logarithmic corrections to the criti- 
cal behavior with the same critical exponents as the pure 
system. Regarding this open question, more recent sim- 
ulations also agree with the exponents being the same as 
the pure model (32j. 



C. Three-dimensional Ising Model 

Figure [7] shows the order-parameter PDF for the three- 
dimensional Ising model universality class according to 
reference [33|. In spite of some progress in obtaining 
this function by analytical methods this is a nice exam- 
ple where simulations are the main source of information 
about the PDF properties. The curve in Figure [7] is given 
by Eq. (3) of reference [33| using the fitting data for the 
lattice 20 3 . 
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FIG. 6: Universal P*(m*) distribution function as function 
of m* = bL^^m for the two-dimensional Ising model with 
spin-1/2, 1, and 3/2 according to Eq. (|34|) . In order to have a 
clearer picture, lesser data than presented in Figure[5]has been 
taken in this case. The extra solid line corresponds to a Wolff 
algorithm to the spin-1/2 model for a lattice size L = 100 
with 10 8 MCS. The error bars are smaller than the symbol 
sizes. 
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the spin-1/2 Baxter- Wu model (when the variables are 
Si = ±1) is depicted in Figure [5] The data in Figure [5] 
have been obtained for a lattice size L — 60 at the exact 
critical temperature fcsT c /J = 2.2692 using Metropolis 
algorithm with 25 million MCS. This distribution is sim- 
ilar to that obtained in reference [35| (the distribution 
shown in Figure [5] has its variance normalized according 
to the present text). Since the Baxter- Wu model corre- 
sponds to three-spin interactions on a triangular lattice, 
the ground state has one ferromagnetic phase with mag- 
netization to = 1 and three ferrimagnetic phases with 
magnetization m = —1/3. This is the reason of having 
one peak at to and another more pronounced peak at 
—m/3. 




FIG. 8: P* distribution function for the two-dimensional 
spin-1/2 Baxter- Wu model. The noise in the data give an 
idea of the error in the simulation. 



bL m 

FIG. 7: Analytical P* distribution function for the three- 
dimensional Ising model according to reference (3^] using the 
results for lattice L = 20. 

D. Spin-1/2 Baxter- Wu Model 

The Baxter- Wu model can be defined by the following 
Hamiltonian 

H = -jJ2SiSjS k , (35) 

where the sum is over all triangles on a triangular lat- 
tice. The exact solution given by Baxter and Wu [34j 
shows a continuous transition at the same critical tem- 
perature as the two-dimensional Ising model with critical 
exponents belonging to the 4-state Potts model univer- 
sality class. The corresponding order-parameter PDF for 



V. EXTENSIVE MIXING FIELDS VARIABLES 

So far we have considered systems presenting symmet- 
ric phase diagrams. In such cases, as for the Ising model 
described by the Hamiltonian (fl~9"f one has the partition 
function 

z= e /3JE < I .J> • s . s j+/ 3ff i:™iS' ! = e tE e +hM e ^ 

states t states £ 

(36) 

where Ei = SiSj and Mi — J2i=i &i are * ne en " 

ergy and magnetization of the labeled state £. We have 
also introduced the variables t = f3J and h = f3H, which 
are the inverse reduced temperature and reduced exter- 
nal field, respectively, where f3 = 1/fceT with ks the 
Boltzmann constant (this definition of /3 should not be 
confused with the order parameter critical exponent pre- 
viously defined). In addition, from above one has the 
mean values 



9 



where the extensive quantities E and M are related to the 
intensive fields t and h, respectively. The corresponding 
free energy is given by G(T, H) = G(t, h) = -k B T\xiZ. 
The second-order phase transition is located at T c and 
H = 0, so for this symmetric problem one has the scaling 
fields (T - T c , H) as in Figured 

In most cases, however, the phase transition line in the 
field space is not symmetric as in the Ising model. For 
instance, consider the Blume-Capel model defined by 



i — sw. 



N 



N 



h = - j E s i s o - H E Si - A E s ^ ( 38 ) 
i=i i=i 

with A being the crystal field and S > 1. In this case 
there is a first-order phase transition line in the temper- 
ature versus crystal field plane (for H = 0) ending at 
a multicritical point |36|. The corresponding phase dia- 
gram regarding the first-order transition looks like that 
sketched in Figure 01 The corresponding partition func- 
tion is given by 



E 



Z = > e r 

states £ 



E 



e 



tE e + hM e +qQi 



states £ 

So, in addition to Eqs. (j3"7| we have here 

1 dZ 



(Q) 



Z dq 



(39) 



(40) 



where q = f3A and (Q) is the mean value of the square 
of the spins, or quadrupole moment, obtained from 



N 



(41) 



Close to the multicritical point (located on the plane 
defined by H = 0) one would take as a possible and naive 
choice the fields (T— T c , A— A c ) where, for simplicity, the 
multicritical coordinates are designated by (T c , A c , H c ). 
This choice, as has been discussed in section [TlJ is not 
convenient and mixing fields have been proven to be more 
desirable to describe the scaling behavior close to the 
continuous transition. 



A. General Approach 

Let us first discuss the general approach to this prob- 
lem which is not applicable only to the Hamiltonian (|38[) 
but to all systems where the first order line is not located 
at constant field A (in this case Q is the corresponding 
extensive quantity of the field q = j3 A). Eqs. (fT"3"|) and 
(|14[) can then be written in a more general form by defin- 
ing t = t — t c and w = q — q c , so that one has 



r) = w + ri, 



(42) 



(43) 



Note that, in this case, close to the multicritical point 
(T c , A c ) instead of the axes T—T c and A— A c as in Figure 
HI we have now t — tt^t- t % t and Gj — which 
differ by a metric factor from the former ones depicted 
in Figure |31 

The question now is how to derive the extensive quan- 
tities associated to the fields ij and r which are not ex- 
plicitly present in the Hamiltonian (f38|) . Nevertheless, 
analogously to Eq. ((3"T]) one can define 



(E) 



1 dZ 
~Z dr' 



1 dZ 



(44) 



and we get 



1 /dZdt 
Z \ dtdr 

W Z \ dt dr, 
Inverting Eqs. (14"2"T) and P3^1 we obtain 

-(rj-rr), 



dZ dti 
dui dr 
dZduj 
dCj drj 



ui 



1 



t 



1 



-(t + sn). 



(45) 



(46) 



(47) 



As the partial derivatives of the partition function with 
respect to t and Cj are identical to the ones with respect 
to i and q, we can use the results (l37l) and (|40"]) to write 
the extensive variables in Eqs. (|45"l) as 



E 



1 



1 + rs 



1 



1 + rs 



(E - rQ), 



(Q + sE). 



(48) 



(49) 



The main task now is to get the probability distribu- 
tion function of these extensive variables. As before, from 
an MC simulation on a finite lattice of linear size L, one 
is able to straightforwardly measure the joint probabil- 
ity distribution Pl(-Ei,Qi) of having discrete values of 
energy Ei and quadrupole moment Qi. From this proba- 
bility distribution we can obtain the corresponding prob- 
ability distribution Pz, (Ej, Qj) for the discrete extensive 
operators Ej and Qj for this lattice size. In fact, what 
we need for this problem is to get the probability dis- 
tribution Vl{Q) of the continuous variable Q obtained 
from the probability distribution of the discrete variable 
Q. The procedure is quite similar to that described in 
the previous section for the magnetization of symmetric 
models. 

In order to do that we first consider a discrete variable 
Xi with normalized probability P(xi) where 



(x) = y j x i P{x i ). 



(50) 
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We also have 



(ax) — aXiP(xi) — a (x) 



(51) 



For the variable yi — axi one has a probability distribu- 
tion P'{yi) in such a way that 



(y) = ax i P'(ax l ) =a(x) . 



(52) 



which, comparing Eqs. (1511) and (1521) . gives P'(axi) = 
P(xi). This means that one is able to get the probability 
distribution Pt(Qi) for the variable Qj defined in Eq. 
(|4*9")) and, as we will see below, this is the key distribution 
for getting the desired universal one. As the MC data 
are stored in pairs of configurations (Ei,Qi), from these 
configurations we can easily construct the probability for 
the composed variables Qi + sEi by just counting their 
corresponding number of occurrences in the total run. 
This probability is equal to the probability of the variable 
Qi. In this way the probabilities Pl(Ei) and Pi(Ej,Qj) 
could also be computed, if necessary. Note still that from 
Eq. (|49|) the variable Qi is discrete, but not integer, 
and not equally spaced, in such a way that we cannot 
use something similar to Eq. (|2"5)) to get its continuous 
distribution. 

On the other hand, we can define the variable 



Qi 



1 {Qj + sE % 

1 ~\- TS ^max 



i 



1 + rs 



(53) 



where Y max is the maximum value of the variable Y^ 
Qi + sEi and yt = (g*±fgi2 . Note further that 



P L (Qi) = PM = P L (Qi), 



(54) 



and yi can be viewed as a continuous variable ranging in 
the interval — 1 < yi < 1 (in general, the mean value of 
Yi is not zero, and we are assuming the distribution is 
dislocated to the right — a different value of Y max can 
be taken when the distribution is dislocated to the left, 
in order to limit the yi variable to the interval —1 to 1). 
Now, for a continuous distribution P(x) one has 



P(x)dx = 1, 



xP(x)dx = (x) . 



(55) 



For the probability distribution V(y) with y — ax one 
has 



V(y)dy 



V(ax)adx = 1, (56) 



which gives 



aV(ax) = P(x). (57) 
From the above equation and Eqs. (l54l and (|53[) one has 
1 



1 + rs 



v l {q) = p L m = p L (Qi)- 



(58) 



We have only worked with the variable of interest, in 
this case Q. One could, however, generalize to the energy 
like variable £ in order to get the distribution P(£,Q) in 
such a way that 

j^—VLiQ) = I Pl(£, Q)d£. (59) 

1 + rs 1 + rs J 

The distribution probability Vl(Q) satisfies the scaling 
relation similar to Eq. (fT8l) . namely 



(60) 



V L {Q)=k Q L^Vl{K Q L^Q), 



where V* is the universal distribution. From above, Eq. 
(l58l) can be written as 



Pt. 



~ -^ M (Q + sE) ), (61) 



1 + rs 



1 + rs 



which is of the same form as Eq. (|2"9"]) . Note that the 
factor 1/(1 + rs) is absorbed in the variables when the 
variance is unit so, analogous to Eq. (f34|) . one gets 



(62) 



P* L ((Q i + sE l )/a))=aP L ( 



where Pi(Qi) is the measured distribution in an MC sim- 
ulation and a is its variance. The temperature, crystal 
field and s are parameters that are tuned in order to 
get the desired universal function. Note, further, that 
(Qi + sEi) is different from zero, so one has also to shift 
the variable to (Qi + sEi) — (Qi + sEi). The distribu- 
tion V* constructed above provides also a method to 
determine the first-order transition line, which happens 
when the shape of V* presents two symmetric peaks at 
the same height. Of course, its shape will not fit any 
universal one when the transition is indeed of first or- 
der, which prevents to misidentifying a continuous tran- 
sition. Besides, analytical continuation of first-order line 
is achieved due to finite-size effects, and critical points 
should be located accordingly, as we will see in the next 
section. 



VI. DISCRETE SPIN AND LENNARD-JONES 
MODELS 

We will present below the use of the mixing fields of 
the previous section to some discrete spin models as well 
as to the liquid-gas critical point. 



A. Two-dimensional Spin-1 Blume-Capel model 

In this case the probability distribution function 
Pl(E,Q) and Pl(Q) can be measured in an MC sim- 
ulation. One can also easily compute Pl(Q + sE) 
and its variance a leading to the desired distribution 
V*((Q + sE)/a)) according to Eq. For a given 

value of the crystal field A/ J one can then tune the tem- 
perature j3J and s in order to have a double peaked 
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function with the same heights. This gives the loca- 
tion of the first-order transition, and from this proce- 
dure the corresponding line, as well as its analytic con- 
tinuation, can be determined. An example of the PDF 
so obtained is given in Figure [9] for the two-dimensional 
spin-1 model, where Si — ±1,0. In order to get the con- 




0.0 

c[(Q+sE)-<Q+sE>] 

FIG. 9: Probability distribution function of the conjugated 
extensive variable for the spin-1 Blume-Capel model (L — 40) 
showing the double peak at the same height. In this example 
the system is, according to reference [37]], at its tricritical 
point. The errors are smaller than the symbol sizes. 

tinuous transition on this first-order line one computes 
the fourth-order cumulant as a function of temper- 
ature and seek for their common crossing. This is ilus- 
tratcd in Figure I1Q[ which gives the tricritical tempera- 
ture and also the tricritical crystal field. The estimates 
are k B T t /J = 0.608(1) and A t /J = 1.9665(3), compara- 
ble to k B T t /J = 0.609(4) and A t /J = 1.965(5) from ref- 
erence ! M and k B T t /J = 0.609(3) and A t /J = 1.966(2) 
from reference [39j. The corresponding order-parameter 
PDF, which is universal for this class of tricritical phe- 
nomena, is shown in Figure [TTJ One can clearly see the 
triple peaks due to the three states of the spin variables. 
To our knowledge no such distribution for the Blumc- 
Capel model in three dimensions has been published in 
the literature. However, the shape for the order parame- 
ter PDF should be qualitatively similar to that shown in 
Figure QT] 



B. Liquid-Gas Critical Point 

The liquid-gas critical point has been analyzed through 
MC simulations by Bruce and Wilding in 1992 [20ll2]j. 
Indeed, it was in these seminal papers that they intro- 
duced the use of the field mixing approach in computer 
simulations for treating asymmetric models and the util- 
ity of the universal probability distribution function given 
by Eq. (j34|) . In references [2(| Hl[ the two-dimensional 
Lennard-Jones fluid was studied through Monte-Carlo 
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FIG. 10: Fourth-order cumulant U% of the Blume-Capel 
model along the first-order transition line and its analytic con- 
tinuation according to reference [371 ]. Errors are not shown; 
only the single histogram result for each lattice is shown. 
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FIG. 11: Probability distribution function of the order pa- 
rameter for the Blume-Capel model at the tricritical point 
for L — 40 according to reference [3?T ]. The errors are smaller 
than the symbol sizes. 



simulations in the grand-canonical ensemble, from which 
the first-order transition line (and its analytical contin- 
uation) has been determined. The dimensionless energy 
$ of the pair interaction has the form 



$ = 4w[((j/r) 12 - (a/rf 



(63) 



where a is a constant that sets the interaction range and 
u! measures the well depth in units of k B T. Differently 
from the previous case of the Blume-Capel model, the 
critical point could be determined in this case by just 
seeking the distribution along the first-order line where a 
match with that of the Ising universality class is achieved. 
The results for this model are 4/w c = 0.440 ± 0.05, /i c = 
-2.20±0.04, and p c = 0.368±0.003. As a result, not only 
the critical point is located but also its universal behav- 
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ior is determined. Figure [T2l shows the PDF at criticality 
compared to that of the spin-1/2 Ising model. This is in 
fact clear evidence, coming from computer simulations, 
of the fluid-magnet universality, which is independent of 
any critical exponent value. In addition, an estimate 
of the ratio of the critical exponents /3/u = 0.125(1) is 
achieved, quite close to the exact value 1/8. A quan- 
titative comparison between the revised scaling method 
proposed by Bruce and Wilding [20] and the complete 
scaling procedure by Kim, Fisher, and Orkoulas [22j can 
be found in reference 12411. 



this probability distribution matches that given in Figure 
[5] Note that here, differently from the tricritical point of 
the spin-1 Blume-Capel model, one does not need to re- 
sort to the cumulant crossings in order to get the contin- 
uous transition. For the present model, instead of having 
a tetracritical point (as should at first sight be expected) , 
one has in fact a double critical end point (two different 
coexisting critical systems), and it was also possible to 
convincingly show that this point indeed belongs to the 
same universality class as the critical ones. 
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FIG. 12: The fluid probability distribution function of the 
density p at the critical point according to references [2(| HJ| . 
The solid line is the order-parameter probability distribution 
of the Ising universality class and the circles and squares are 
the data for the Lennard- Jones fluid for lattice sizes L = 16 
and L = 8, respectively. Acj is the energy deviation. 

Three-dimensional fluid models have also been ana- 
lyzed according to the present approach. The extension 
to the three-dimensional Lennard- Jones fluid has been 
treated by Wilding [40j, and the decorated lattice gas 
and a polymer system have been studied by Wilding and 
Miiller [4J|. In all cases, the data collapsed onto the 
three-dimensional PDF shown in Figure [7J providing the 
location of the critical point as well as the corresponding 
expected universality class. 



C. Spin-3/2 Blume-Capel model 



The Blume-Capel model Hamiltonian (155)) for spin-3/2 
means that Si — ±3/2, ±1/2. Besides a second-order 
phase transition line, the corresponding phase diagram 
shows a first-order line ending up at a double critical end 
point [42|]. Within the present approach, the first-order 
line (which is in fact a quadruple line where four phases 
are coexisting) can be determined when the extensive 
mixing field variables' probability distribution function 
has peaks of equal heights, and the double critical end 
point can be located, with a reasonable precision, when 



D. Spin-1 Baxter- Wu model 

The spin-1 Baxter- Wu model is a generalization of the 
Hamiltonian (|35[) by considering S% = ±1,0 and with a 
crystal field interaction in such a way that the Hamilto- 
nian reads 



A' 



U = - J^2 SiSjSk - A ^ Sf 



(64) 



i,j,k 



The phase diagram of the above model is similar in shape 
to the spin-1 Blume-Capel model. However, instead of 
a triple line ending at a tricritical point one expects a 
quintuple line (five phases coexisting) ending at a pen- 
tacritical point [43]. The second-order line should be in 
the same universality class as the spin-1/2 case and the 
pentacritical point belongs to a different class of univer- 
sality. Figure [13] depicts the order-parameter PDF for 
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FIG. 13: P* distribution function of the spin-1 Baxter- Wu 
model. The values of A/ J correspond to the continuous tran- 
sition and are taken at reduced temperatures ksT/J where 
the distribution has the same form as that of the spin-1/2 
model from Figure [8] (which is also plotted in this Figure) . 
The inset shows the corresponding distribution function at 
the pentacritical point. The noise in the distributions give an 
idea of the errors from the simulations. 

values of A/ J in the second-order line together with the 
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distribution of the spin- 1/2 case [Hj]. One can see that 
the continuous transition is in fact in the same universal- 
ity class, as expected. The inset shows the corresponding 
distribution function at the multicritical point. As in this 
case we also have a paramagnetic phase the distribution 
shows an extra peak for small values of the magnetization 
M. 



VII. SOME RELATED MODELS AND FINAL 
REMARKS 

Other related models have also been studied according 
to the present approach. For instance, regarding discrete 
state systems, the isolated critical point of the asymmet- 
ric Ising model with two- and three-spin interactions on a 
triangular lattice, in the presence of an external magnetic 
field, has been obtained and it has been shown to belong 
to the Ising universality class [45|,|46|. On the other hand, 
the spin- 1/2 Ising model on the square lattice with two- 
body and three-body interactions has been studied [47]] 
and very precise distribution functions for much larger 
systems show the evidence for the "universality" of the 
scaled distribution, including a tricritical point. It is in- 
teresting, in this particular case, to make a comparison of 
the probability distribution function of the order param- 
eter shown in Figure 111! with the corresponding one for 
the tricritical point in the Ising model given in Figure 6 
of reference [47| • While in the former case we have a cen- 
tral peak which is more intense than the two side peaks 
(because both the Si = and disordered states have 
zero magnetization), for the tricritical point in the Ising 
model, all three peaks have the same intensity (because 
there is no Si = state). Still regarding discrete mod- 
els, the Q-state Potts model has also been investigated 
through this technique. The large-Q Potts model in an 
external field presents an isolated critical point which 
can also be located through the use of the mixing fields 

[Mill- 
Continuous spin systems and </> 4 models have also been 
treated by the mixing fields procedure. The PDF of the 
XY model in two dimensions has been measured [49[ as 
well as for the isotropic Heisenberg model. Isolated crit- 
ical points can be achieved when considering the vector- 
ized Blume-Emery-Griffiths model applied to the phase 



transition in He 3 -He 4 mixtures [13, HJ. The PDF of the 
(j) 4 model has also been studied in reference 27]. 

A different algorithm, namely the grand canoni- 
cal Monte Carlo method, has been employed with 
the present scaling technique to study confined lattice 
homopolymers[52j. In a quite recent paper (53[, the field 
mixing revised scaling has been used in order to study 
the fluid phase behavior of an athermal model of colloids 
and nano absorbing polymers on a finite lattice. More re- 
cently, an extension of the probability distribution func- 
tion for systems in non equilibrium has also been devised 
in the case of the contact process problem [54| . 

Thus, it seems that the revised scaling proposed by 
Rehr and Mermin, associated with the mixing field recipe 
developed by Wilding and Bruce, is still a useful tool 
in treating asymmetric magnetic models exhibiting first- 
order lines and multicritical points. Naturally, a com- 
parison study by applying the complete scaling proposed 
by Kim and Fisher on magnetic systems would be very 
welcome. 
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